Obiektem analizy był zbiór danych dotyczących połowu śledzia oceanicznego. Dane zwierały ponad 50 tysięcy wpisów zebranych na przestrzenii około 60 lat. Podczas analizy zbioru zauważono, że od pewnego momentu długość śledzi zaczęła regularnie spadać. Po wstępnym oczyszczeniu danych i usunięciu z nich trywialnych zależności, wykorzystano algorytm uczenia maszynowego eXtreme Gradient Boosting i wykryto, że główną przyczyną tego zjawiska jest utrzymujący się wzrost temperatury przy powierzchni wody.
library(dplyr)
library(ggplot2)
library(grid)
library(gridExtra)
library(reshape2)
library(caret)
library(knitr)
library(plotly)
set.seed(666)
Dane zostały wczytane z pliku CSV. Brakujące wartości oznaczone były znakiem ?, który zamienono na wartość NA.
data <- read.csv("sledzie.csv")
data[data=="?"] <- NA
data <- mutate_all(data, as.numeric)
completeData <- data[complete.cases(data),]
Zbiór danych zawierał 52582 wpisów, każdy z nich opisany był przez 16 atrybutów numerycznych, w tym jedną kolumnę identyfikatora, porządkującą zbiór chronologicznie.
Krótkie podsumowanie najważniejszych statystyk dla każdego z atrybutów (wyłączono atrybut identyfikatora X).
knitr::kable(summary(data[,2:9]))
| length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | fbar | |
|---|---|---|---|---|---|---|---|---|
| Min. :19.0 | Min. : 2.00 | Min. : 2.00 | Min. : 2.00 | Min. : 2.00 | Min. : 2.00 | Min. : 2.00 | Min. :0.0680 | |
| 1st Qu.:24.0 | 1st Qu.: 2.00 | 1st Qu.:14.00 | 1st Qu.:13.00 | 1st Qu.:15.00 | 1st Qu.:13.00 | 1st Qu.:13.00 | 1st Qu.:0.2270 | |
| Median :25.5 | Median :15.00 | Median :24.00 | Median :20.00 | Median :27.00 | Median :23.00 | Median :27.00 | Median :0.3320 | |
| Mean :25.3 | Mean :16.39 | Mean :23.64 | Mean :23.45 | Mean :28.07 | Mean :23.85 | Mean :27.82 | Mean :0.3304 | |
| 3rd Qu.:26.5 | 3rd Qu.:28.00 | 3rd Qu.:34.00 | 3rd Qu.:36.00 | 3rd Qu.:41.00 | 3rd Qu.:35.00 | 3rd Qu.:43.00 | 3rd Qu.:0.4560 | |
| Max. :32.5 | Max. :40.00 | Max. :49.00 | Max. :49.00 | Max. :52.00 | Max. :49.00 | Max. :52.00 | Max. :0.8490 | |
| NA | NA’s :1581 | NA’s :1536 | NA’s :1555 | NA’s :1556 | NA’s :1653 | NA’s :1591 | NA |
knitr::kable(summary(data[,10:16]))
| recr | cumf | totaln | sst | sal | xmonth | nao | |
|---|---|---|---|---|---|---|---|
| Min. : 140515 | Min. :0.06833 | Min. : 144137 | Min. : 2.00 | Min. :35.40 | Min. : 1.000 | Min. :-4.89000 | |
| 1st Qu.: 360061 | 1st Qu.:0.14809 | 1st Qu.: 306068 | 1st Qu.:14.00 | 1st Qu.:35.51 | 1st Qu.: 5.000 | 1st Qu.:-1.89000 | |
| Median : 421391 | Median :0.23191 | Median : 539558 | Median :25.00 | Median :35.51 | Median : 8.000 | Median : 0.20000 | |
| Mean : 520367 | Mean :0.22981 | Mean : 514973 | Mean :25.45 | Mean :35.51 | Mean : 7.258 | Mean :-0.09236 | |
| 3rd Qu.: 724151 | 3rd Qu.:0.29803 | 3rd Qu.: 730351 | 3rd Qu.:36.00 | 3rd Qu.:35.52 | 3rd Qu.: 9.000 | 3rd Qu.: 1.63000 | |
| Max. :1565890 | Max. :0.39801 | Max. :1015595 | Max. :52.00 | Max. :35.61 | Max. :12.000 | Max. : 5.08000 | |
| NA | NA | NA | NA’s :1584 | NA | NA | NA |
Pomimo dużej liczby danych w zbiorze, ich zróżnicowanie było bardzo niewielkie. Wynikało to z faktu, że w większości atrybuty opisują wartości mierzone w oknie rocznym (dla danego roku parametry się nie zmieniają). Brakujące wartości w odpowiednich kolumnach zostały zastąpione wartościami z odpowiednich kolumn pełnych rekordów, w których wartość atrybutu recr była taka sama (atrybut ten opisuje roczny narybek, więc dla danego roku jest on stały).
W celu przyśpieszenia wyszukiwania odpowiednich wartości zastępujących, dane zostały pogrupowane wzlędem kolumny recr, a następnie umieszczone w tablicy haszowej. Przedstawiona poniżej funkcja replace_missing zawiera instrukcję warunkową, dla przypadku, gdy dany rekord nie ma atrybutu length - ten przypadek jest zastępowany średnią długością z danego roku i miesiąca. W dostarczonym ziobrze danych, taki przypadek jednak nie wystepował.
replacementData <- group_by(completeData, recr) %>% filter(row_number()==1) %>% ungroup() %>% arrange(recr)
lookupReplacement <- new.env()
for(ridx in 1:nrow(replacementData)){
key <- as.character(unlist(replacementData[ridx,"recr"]))
lookupReplacement[[key]] <- replacementData[ridx,]
}
replace_missing <- function(dataRow) {
recrid <- unlist(dataRow["recr"])
replacement <- lookupReplacement[[as.character(recrid)]]
dataRow[which(is.na(dataRow))] <- replacement[which(is.na(dataRow))]
if(is.na(dataRow["length"])) {
m <- unlist(dataRow["xmonth"])
l <- unlist((filter(completeData,recr==recrid,xmonth==m) %>% group_by(recr) %>% summarize(l=mean(length)))["l"])
dataRow["length"] <- l
}
dataRow
}
withMissing <- data[!complete.cases(data),]
if(nrow(withMissing) > 0){
for(item in 1:nrow(withMissing)){
withMissing[item,] <- replace_missing(withMissing[item,])
}
}
Przed rozpoczęciem analizy statystycznej, dane we wszystkich kolumnach (z wyjątkiem kolumny identyfikatora X) zostały znormalizowane do przedziału od 0.0 do 1.0.
completeData <- arrange(merge(completeData,withMissing,all=TRUE),X) %>% mutate_all(as.numeric)
normalized <- completeData
for(col in colnames(normalized)){
if(col == "X")next
normalized <- mutate_(normalized,.dots=setNames(paste0("(",col,"-min(",col,"))/(max(",col,")-min(",col,"))"),col))
}
Korelacja pomiędzy atrybutami została sprawdzona za pomocą współczynnika korelacji Pearsona. Wykres przedstawia zależności pomiędzy atrybutami (ciemniejszy kolor, oznacza silniejszą korelację). Dla ułatwienia interpretacji, wykres przedstawia wartości bezwzlędne współczynnika korelacji.
Najsilniesze powiazanie występowało pomiędzy atrybutami:
Poniższy animowany wykres przedstawia jak zmieniał się trend długości śledzi z ubiegiem czasu - od początku do końca pomiarów. Łatwo zauważyć, po około 2/5 czasu uwzlędnionego w zbiorze nastąpiło przełamanie trendu i długość śledzi zaczeła maleć.
Wykres został wygenerowany za pomocą kodu w komentarzu, RMarkdown miał problem z użyciem bilbioteki animate i wykres nie generował się poprawnie
Dla łatwiejszej analizy, poniższe wykresy przedstawiają trend dla kilku kluczowych punktów:
library(animation)
show_on<-c(20,40,60,100)
max_x <- max(completeData[,"X"])
x_es <- 1:max(max_x)
x_es <- x_es[x_es%%(floor(length(x_es)/100))==0]
trend_plots <- c()
for(pseudo_x in x_es){
percent <- ceiling((pseudo_x*100)/max_x)
x_label <- paste0("czas, który upłynął: ",percent," %")
p <- ggplot(filter(completeData, X<=pseudo_x),aes(X,length)) +
geom_smooth(method="gam",formula = y ~ s(x, bs = "cs")) + labs(y="trend długości śledzi [cm]",x=x_label) +
ylim(25,28) +
theme( axis.ticks.x = element_blank(), axis.text.x = element_blank())
trend_plots <- c(trend_plots, p)
if(any(show_on==percent)){
print(p)
}
}
saveGIF({
for(idx in seq_along(trend_plots)){
print(trend_plots[idx])
}
}, interval = 0.2, ani.width = 550, ani.height = 350)
Poniższe histogramy prezentują rozkłady wartości atrybutów w zbiorze. Poza atrybutem length, który charakteryzuje się rozkładem zbliżonym do rozkładu normalnego, atrybuty posiadają bardzo nieregularne rozkłady wartości, z kilkoma wartościami odstającymi. Nie pokazano rozkładu dla atrybutów X (identyfikator) oraz xmonth.
Tak jak wspomniano wcześniej (podpunkt Brakujące Wartości) większość atrybutów zawiera mało unikalnych wartości. Warto zwrócić uwagę na oś Y - liczba elementów, której skala rozciąga się od 0 do 6000. Dla niektórych wartości poszczególnych atrybutów, ich powtarzalność sięga prawie maksymalnej wartości tego przedziału.
Przed rozpoczęciem trenowania modelu, dane zostały znormalizowane do przedziału od 0.0 do 1.0 (wszystkie atrybuty, z wyjątkiem długości - length, który zachował oryginalne wartości). Usunięto również atrybuty o silnej korelacji (na podstawie analizy opisanej wcześniej). Usunięte atrybuty to: chel1, chel2, cumf oraz X (kolumna identyfikatora).
Zbiór został podzielony na treningowy i testujący w stosunku odpowiednio 7:3, przy czym zastosowano losowanie warstowe, zachowujące rozkład wartości atrybutu length. Po wielu eksperymentach, najlepiej ocenianym modelem predykcyjnym okazał się eXtreme Gradient Boosting. Jego parametry zostały dobrane za pomocą metody kroswalidacji.
withoutCorelated <- select(normalized,-chel1,-chel2,-cumf, -X)
herrings <- withoutCorelated
herrings[["length"]] <- completeData[["length"]]
inTraining <- createDataPartition(herrings$length, p=0.7, list=FALSE)
training <- herrings[inTraining,]
testing <- herrings[-inTraining,]
ctrl <- trainControl(method="repeatedcv",number=3,repeats=3)
model <- train(length ~ .,data = training, method="xgbLinear", trControl=ctrl)
Najlepsze parametry oraz miary oceny dla zbioru treningowego:
b<-model$bestTune
kable(b)
| nrounds | lambda | alpha | eta | |
|---|---|---|---|---|
| 19 | 50 | 0.1 | 0 | 0.3 |
best<-filter(model$results,nrounds==b$nrounds&lambda==b$lambda&alpha==b$alpha&eta==b$eta)
kable(select(best,RMSE,Rsquared))
| RMSE | Rsquared |
|---|---|
| 1.143814 | 0.5228055 |
Na danych testowych model uzyskał następujące miary:
predictions <- testing
predictions[["predicted"]] <- predict(model, newdata=testing)
rmse_testing <- sqrt(mean((predictions$length-predictions$predicted)^2))
ss_res<-sum((predictions$length-predictions$predicted)^2)
ss_tot<-sum((predictions$length-mean(predictions$length))^2)
rsq_testing <- 1 - (ss_res/ss_tot)
kable(data.frame("RMSE"=rmse_testing,"Rsquared"=rsq_testing))
| RMSE | Rsquared |
|---|---|
| 1.131134 | 0.5280382 |
p <- ggplot(predictions,aes(length,predicted)) + geom_point() + geom_smooth(method="glm") + labs(y="długość przewidziana przez model [cm]",x="rzeczywista długość [cm]") + xlim(19,32.5) + ylim(19,32.5) + theme_minimal()
ggplotly(p)
Ranking ważności atrybutów dla zbudowanego modelu regresji przedstawia się następująco:
importance <- varImp(model)
ggplot(importance,aes(Overall)) + geom_bar(stat="identity",fill="#0089B2") + theme_minimal() + labs(y="ważność", x="atrybut",title="Ranking atrybutów")
Najważniejszym atrybutem okazal się sst, oznacza to, że głownym czynnikiem wpływającym na długość śledzia jest temperatura przy powierzchni wody. W początkowym okresie, gdy temperatura spadała, średnia długość śledzi wzrastała. Po osiągnięciu minimum temperatury, zaczęła ona dość gwałtownie rosnąć i od tego samego punktu, długość śledzia zaczęła nieustannie maleć. Na poniższym wykresie można zauwazyć nietrywialną korelację pomiędzy tymi atrybutami:
ggplot(normalized,aes(X,sst)) + geom_smooth(aes(X,sst,colour="temperatura (sst)")) + geom_smooth(mapping=aes(X,length,colour="długość śledzia")) + labs(title="Korelacja sst i length", y="znormalizowana wartość",x="czas") + theme(legend.text =element_text("atrybut")) + scale_colour_manual("",values=c("#0089B2","red")) + theme_minimal() + theme(legend.position = "bottom")
Pozostałe atrybuty o wiele mniej wpływały na długość śledzi, jednak na niektórych z nich (np. dostępnośc planktonu) można również zauważyć korelację ich trendu względem długości śledzi.
Po odkryciu wspomnianych zależności, usunięto 4 najważniejsze atrybuty ze zbioru (sst, recr, totaln, lcop1), aby sprawdzić, jak zachowa się model. Okazało się, że zbudowany model ma bardzo zbliżoną jakość, pomimo wydawałoby się “trudniejszego” zbioru danych. Tym razem, ważność atrybutów jest nieco bardziej zbilansowana - poza najistotniejszym - natężeniu połowów (fbar), atrybuty sal, nao oraz cfin1 są na podobnym poziomie istotności.
herrings2 <- select(withoutCorelated,-lcop1, -sst, -totaln, -recr)
herrings2[["length"]] <- completeData[["length"]]
inTraining2 <- createDataPartition(herrings2$length,p=0.7,list=FALSE)
training2 <- herrings2[inTraining2,]
testing2 <- herrings2[-inTraining2,]
ctrl <- trainControl(method="repeatedcv",number=3,repeats=3)
model2 <- train(length ~ .,data = training2, method="xgbLinear", trControl=ctrl)
predictions2 <- testing2
predictions2[["predicted"]] <- predict(model2, newdata=testing2)
rmse_testing2 <- sqrt(mean((predictions2$length-predictions2$predicted)^2))
ss_res2<-sum((predictions2$length-predictions2$predicted)^2)
ss_tot2<-sum((predictions2$length-mean(predictions2$length))^2)
rsq_testing2 <- 1 - (ss_res2/ss_tot2)
importance2 <- varImp(model2)
ggplot(importance2,aes(Overall)) + geom_bar(stat="identity",fill="#0089B2") + theme_minimal() + labs(x="atrybut", y="ważność",title="Ranking atrybutów modelu alternatywnego")
Wskaźniki dla alternatywnego modelu dla danych treningowych:
| RMSE | Rsquared |
|---|---|
| 1.146805 | 0.5211433 |
oraz testowych:
| RMSE | Rsquared |
|---|---|
| 1.143814 | 0.5228055 |
Wejściowy zbiór danych okazał się dość trudny w przetwarzaniu. Dla wielu różnych długości śledzia, wszystkie pozostałe atrybuty posiadały te same wartości, co przyczyniło się do dużych odchyleń w wartościach przewidywanych przez budowane modele predykcyjne.
W trakcie analizy, zgodnie z podanymi informacjami, traktowano X jako atrybut wskazujący na czas, pomimo braku jednostki i jego niespójności ze zmianami wartości atrybutu xmonth. Przy analizie tego typu danych wejściowych, warto poprosić dostawcę danych o przekazanie dodatkowych atrybutów, jeżeli jest taka możliwość.
Wygenerowanie całego raportu od podstaw zajmuje około 5 minut na procesorze i5 @ 1.8GHz